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This paper presents the results of a one-dimensional numerical simulation of the transient chilldown of a 
vertical stainless steel tube with liquid nitrogen. The direction of flow is downward (with gravity) through the 
tube. Heat transfer correlations for film, transition, and nucleate boiling, as well as critical heat flux, 
rewetting temperature, and the temperature at the onset of nucleate boiling were used to model the 
convection to the tube wall. Chilldown curves from the simulations were compared with data from 55 recent 
liquid nitrogen chilldown experiments. With these new correlations the simulation is able to predict the time 
to rewetting temperature and time to onset of nucleate boiling to within 25% for mass fluxes ranging from 
61.2 to 1150 kg/m 2 s, inlet pressures from 175 to 817 kPa, and subcooled inlet temperatures from 0 to 14 K 

below the saturation temperature. 
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Nomenclature 

= cross-sectional area of tube material 
= specific heat at constant pressure 
= inner diameter of tube 
= two-phase Reynolds number factor 
= mass flux 

= heat transfer coefficient 
= thermal conductivity 
= Nusselt number, hD/k 
= pressure 

= Prandtl number, C r ,p!k 
= heat rate 
= heat flux 

= Reynolds number, GD/p 
= suppression factor 
= time 

= tube thickness 
= temperature 

= temperature at onset of nucleate boiling 
= wall temperature 
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time step 
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convection to the wall from inside fluid 

critical heat flux 

from data measurements 

with respect to inner diameter 

Dittus-Boelter 

fluid 

film boiling 
node number 

liquid property at saturation temperature 

with respect to length downstream of inlet 

nucleate boiling 

parasitic heat 

solid tube wall 

saturation property 

from the simulation 

transition boiling 

vapor property at saturation temperature 

wall 


I. Introduction 

The use of cryogenic propellants in chemical rocket engines like liquid oxygen (LOX) and liquid hydrogen 
(LH 2 ) is desirable because cryogens are non-toxic and have a higher specific impulse compared to other propellants. 
LH 2 is also used to power nuclear rocket engines which have a significantly higher specific impulse than any 
chemical propellant combination. On the other hand, the transfer of a liquid cryogen is much more difficult than the 
transfer of storable propellants (propellants that maintain a liquid phase at normal temperatures like RP-1 and 
hypergols) in terms of both safety and efficiency. The transfer of storable propellants from a tank to an engine or 
receiver tank is trivial in comparison to the transfer of cryogenic propellants. Cryogens are stored at extremely low 
temperatures, and when they are introduced to a normal temperature transfer line they immediately absorb heat from 
the tube and proceed to boil. This is problematic in most cases since the transfer of vapor-free propellant is critical, 
like in the case of transferring propellant from a propellant storage tank to another tank, engine, or a fuel depot. 
Some quantity of the cryogen must be used to chill down the tube before single-phase liquid can be transferred 
through the entire line. The amount that is used is dependent on the flow rate and heat flux from the tube wall to the 
fluid. 

Unfortunately, the robust, accurate prediction of two-phase flow boiling heat transfer continues to be a 
challenge. The diversity of two-phase flow features that occur during chilldown makes it a complex task to provide 
correlations that are valid over a wide range of conditions. Typically, separate heat transfer correlations are used for 
different boiling regimes and sometimes different flow patterns. It is vital for designers of cryogenic propellant feed 
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systems both on the ground and in microgravity that robust heat transfer correlations are provided to enable accurate 
predictions of the chilldown process. 

This paper presents an experimentally validated one-dimensional (ID) numerical simulation of the chilldown of 
a vertically-aligned tube with liquid nitrogen (LN 2 ) using a new set of two-phase flow boiling heat transfer 
correlations. Brief summaries are first given on the chilldown process with a focus on the different boiling regimes, 
previous cryogenic chilldown models, and the experiment from which the data is used to validate the heat transfer 
correlations used in the current work. Then the heat transfer correlations used in the simulation are described for 
each boiling regime and a description is given on the setup of the numerical simulation. Lastly, the results of the 
simulation are compared to the data. The success of this simulation demonstrates that these correlations can be 
integrated into any one-dimensional, one-fluid model to enable the accurate and economical prediction of transfer 
line chilldown. 


II. Cryogenic Chilldown Process 

The cryogenic chilldown process can be perceived as a reverse boiling curve, which divides the boiling process 
into three different regimes. 1 ' 2 Data from a LN 2 chilldown test at an average mass flux of 39 kg/m 2 s and average 
pressure of 202 kPa is shown in Figure 1, which consists of a temperature vs. heat flux plot on the left and a 
temperature vs. time plot on the right. At the start of the transfer, the tube is at such a higher temperature (between 
270 K and 300 K depending on the ambient conditions) than the saturation temperature of the cryogen (between 77 
K and 110 K for LN 2 depending on the fluid pressure) that any liquid approaching the wall is vaporized, leaving 
only vapor touching the wall. This is called film boiling and it persists for a majority of the chilldown time. Only 
once the temperature of the wall drops below the rewetting temperature (about 120 K for the test in Figure 1) does a 
significant amount of liquid begin to touch the wall. This rewetting temperature is called the Leidenfrost temperature 
for droplet vaporization experiments. This marks the initiation of transition boiling, which consists of unpredictable, 
intermittent intervals of liquid-to-wall contact and vapor-to-wall contact. Transition boiling can be seen as the 
mixture of film boiling and nucleate 
boiling. The heat flux increases as the 
temperature drops further below the 
rewetting temperature until the critical 
heat flux (CHF) is reached, which marks 
the start of full nucleate boiling. During 
nucleate boiling, liquid remains in 
contact with the wall while vapor 
bubbles or slugs still subsist within the 
flow. The CHF is usually the largest 
amount of heat flux during the entire 
chilldown process, and as the 
temperature of the tube wall continues to 
decline, the heat flux decreases until all 
bubbles are condensed. If the 
temperature of the incoming liquid is 
below the temperature at the onset of 
nucleate boiling (ONB), there is not 
enough superheat to sustain nucleate 
boiling. Single-phase convection heat transfer between the wall and the liquid is the only mode of heat transfer left. 
Typically the ONB occurs between 5 K and 10 K above the saturation temperature for LN 2 . 3 

Prediction of the heat flux during the chilldown process is a difficult task because it depends on the mass flux, 
the fluid properties, the tube properties, the amount of gravity, the direction of flow with respect to gravity, and even 
the surface finish of the tube, among other things. Additionally, chilldown is a transient process with a degree of 
unpredictability due to the instabilities of two phase flows and the complicated interaction at the solid-liquid-vapor 
interface. To achieve the best accuracy, heat transfer coefficient (HTC) correlations are provided separately for each 
boiling regime to predict the heat transfer from the wall to the fluid. Moreover, models for the rewetting 
temperature, the CHF, and the ONB temperature are required to determine the dividing points between regimes. 



flux plot on the left, temperature vs. time plot on the right. The different boiling 
regimes are shown. 
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III. Other Chilldown Models 

Many previous simulations of the cryogenic chilldown of tubes have only used one HTC correlation for all three 
boiling regimes. For example, Cross et al., 4 Majumdar and Steadman, 5 and Majumdar and Ravindran 6 each used the 
modified Miropolskii correlation for all two-phase flow conditions. 7 This worked well for predicting chilldown rates 
with hydrogen at high flow rates in long transfer lines (about 60 m). The correlation, however, does not work well 
for other fluids like LN 2 where there is a clear presence of three distinct boiling regimes. Nor does it work well for 
situations of nonequilibrium and flow instability, which occur frequently in cryogenic flow systems. Nonequilibrium 
in flow boiling typically occurs in two situations. One is when subcooled liquid is introduced into a tube at relatively 
high mass fluxes in which the tube temperature is higher than the rewetting temperature. In this situation a vapor 
annulus at the saturation temperature surrounds a subcooled liquid core. Because the two phases of the fluid are at 
different temperatures the system is said to be in thermal nonequilibrium. Calculating the equilibrium quality based 
on the enthalpy can lead to a negative value. The second situation of nonequilibrium occurs in high quality flow in 
which the vapor becomes superheated even with the presence of liquid droplets that are at the saturation 
temperature. Calculating the equilibrium quality based on the enthalpy can lead to a value greater than one. 

Agrawal et al. N simulated the chilldown of a tube with LN 2 using a separate correlation for film boiling and 
nucleate boiling. The Miropolskii correlation was used for film boiling and the Chen correlation was used for 
nucleate boiling. Again, nonequilibrium situations could not be modeled by using the Miropolskii correlation. It was 
not specified what was used for transition boiling. The predicted wall temperatures were not compared to data. Y uan 
et al. 9 developed a two-fluid model for film boiling and used the Chen correlation for both transition and nucleate 
boiling. The two-fluid model required solving two sets of conservation equations for the fluid with several semi- 
empirical closure models. The model worked well against data and took into account nonequilibrium, but required 
careful attention to the numerical solution procedure and very small time stepping to ensure stability and accuracy. 
Accurate correlations that can be applied to a more efficient one-fluid model, but do not sacrifice the ability of 
predicting the heat flux in nonequilibrium conditions, like the correlations of the current work, would be preferred 
over this approach for a ID, one-fluid model simulation. The correlation set presented in this paper goes further than 
other one-fluid cryogenic chilldown models by: 

• applying a new film boiling heat transfer correlation that is usable in nonequilibrium conditions 

• applying a separate correlation for transition boiling 

• giving explicit correlations for CFIF, rewetting temperature, and the ONB temperature to allow easy coding 
logic to determine which FITC correlation to use 


IV. Chilldown Experiment Overview 

LN 2 chilldown experiments were designed, built, and performed at the University of Florida Department of 
Mechanical and Aerospace Engineering in Gainesville, FL. The purpose of the experiment was to obtain chilldown 
curves on a vertically aligned stainless steel tube with a downward flow direction over a wide range of mass flow 
rates and pressure conditions in order to support the development of robust heat transfer correlations in each of the 
three boiling regimes. LN 2 was selected as the test fluid in the present study not only for the safety consideration but 
also because its properties such as the surface tension, boiling point and evaporative latent heat are similar to those 
of LOX. The initial temperature of the test section was kept at normal room temperature, near 293 K, for all the 
runs. 

The tests were operated as follows: At the start, the hardware leading up to the test section was pre-chilled by 
flowing LN 2 from a pressurized gas cylinder from the cryogenic Dewar through a subcooling pool cooler and 
eventually to vent. The subcooling pool cooler gave a subcooling margin between 0 and 14 K at the bypass valve 
just before the test section. After the line upstream of the test section was completely chilled down, which was 
verified by a TC embedded in the flow just upstream of the bypass valve, the three-way bypass valve was opened to 
allow flow into the test section. The test section was housed inside a stainless steel chamber which provided a 
vacuum environment of 0.014 Pa pressure to insulate the test section. The transient temperature and pressure data 
along the tube were collected by T-type thermocouples (TCs) and pressure transducers (PTs). Downstream of the 
test section, the two-phase flow entered a heat exchanger, where the two-phase flow fully evaporated into warm 
nitrogen vapor, which enabled a single-phase flow rate measurement with a gas flow rate meter. The test was 
terminated once the temperature readings on the test section showed that the temperature had dropped well below 
the ONB temperature and maintained a steady condition. 
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Figure 2 shows a schematic of the test section (rotated 90° counterclockwise), which is a 57.2 cm long, 1.27 cm 
outside diameter, 0.0508 cm thick stainless steel tube. For each test, single-phase LN 2 was supplied to the test 
section inlet (just downstream of the bypass valve), at a certain flow rate and pressure. The vacuum chamber was 
used to remove most of the surrounding air in the test section and thus significantly reduced the amount of heat 
transfer from the ambient environment to the tube. This lessened the uncertainty in determining the convective heat 
flux between the inner wall of the test section and the fluid. A set of three TCs was placed at two axial locations on 
the test section inside the vacuum chamber. TC-1, 2, and 3 were located upstream of TC-4, 5, and 6. The TCs were 
attached to the outer tube wall and placed 90° apart around the tube circumference. Two PTs, PT-1 and PT-2, were 
also placed inside the test section 3 cm apart from the TCs. 


Vacuum Chamber 



Figure 2. Test section geometry. TC-1, TC-2, and TC-3 are the three thermocouples that make up the 
upstream TC station. TC4-6 make up the downstream TC station. PT-1 and PT-2 are the upstream and 
downstream pressure transducers, respectively. The bold black box represents the vacuum chamber. 


V. Heat Transfer Correlations 

In this section each FITC correlation for the three boiling regimes is presented, along with correlations for the 
rewetting temperature, the CHF, and the ONB temperature. These correlations are used in the chilldown simulation 
presented in the following sections. 

A. Film Boiling HTC Correlation 

A major portion of LN 2 chilldown is spent in fdm boiling. 12 During film boiling, the flow features can change 
depending on the mass flux, the distance from the tube inlet, etc. For an intermediate to high mass flux a liquid core 
can penetrate from the quenching front to a long distance downstream. The vapor forms an annulus around the liquid 
so that there is no liquid-to-wall contact. This is called inverted annular flow. For low Re numbers this core can be 
negligibly small or nonexistent and the fdm boiling regime consists of a high quality mixture where the liquid is in 
the form of droplets or small slugs. This type of flow is called dispersed flow. These are the two regimes of flow in 
fdm boiling. Several studies have divided these flow regimes into smaller subregimes. This paper will, however, 
utilize a correlation that is employed for the entire film boiling regime, regardless of the flow structure. 

The correlation for the fdm boiling Nusselt (Nu) number, fit to all of the vertical downward flow data, is given 
as: 


Nit/,, = h jh DI k v = (7.55 xlCT 4 — 7.43 xlO^z/ Z))ite v ° 941 (l — x e ) 523 Pr° A + 0.0568 (£, / k r )We z 0j b (1) 

where We D = G 2 D/(p/a) is based on tube diameter D, and dp, = (300 - 7j,)/(300 - T wc j) is the nondimensional fdm 
boiling temperature. The temperatures are in units of Kelvin. The subscript “v” and denote that the property is a 
vapor property and liquid property, respectively, evaluated at the saturation temperature. 

The first of the two terms in this correlation accounts for the heat transfer between the wall and the vapor. It is 
very similar to the single-phase Dittus-Boelter equation except that it is dependent on the distance from the inlet, and 
that it includes two-phase effects. The farther from the inlet, the smaller the heat transfer to the vapor which is likely 
a result of the developing thermal and hydrodynamic boundary layers near the inlet. The (1 - x e ) term accounts for 
the size of the vapor layer and thus the actual speed of the vapor layer. As the quality increases the vapor layer 
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thickens, and from continuity, the velocity of the vapor near the wall increases. The term (1 —x e ) was chosen instead 
of including x e inside Re v because in the case of a subcooled inlet where x e is negative; (1 - x e )~ 5 23 produces a real 
number, but including a negative x inside Re v would produce an imaginary number when taken to the power of 
0.941. Pr v accounts for the differences between fluids. The exponent is left as 0.4 to match the single-phase Dittus- 
Boelter equation. 

The second term in Eq. (1) represents the added heat transfer from liquid impingement against the wall. The 
breakup of the liquid core and subsequent two-phase mixing can be so violent that it causes liquid slugs or droplets 
to eject from the liquid core and impinge on the wall and vaporize, effectively increasing the heat flux from the wall. 
This effect is more significant for higher mass fluxes. It has been shown that the breakup of the liquid core is a 
strong function of the We number. 10 ' 12 From the current work it was estimated that the added heat transfer from this 
effect increases linearly with We number. It was also determined from analyzing the data that the liquid 
impingement effect increases nonlinearly with 9^. 0 /h goes to zero at T w = 300 K. (The constant in 0 /h was chosen as 
300 K instead of the test starting temperature of 293 K because it was a round number to base the correlation on. In 
this case when T w is 293K 9 fi is still almost zero and can be considered negligible.) For T w > 300 K, 6 ^ should be set 
to zero, although it would be wise for the user to validate the correlation at temperatures above this temperature. As 
T w is lowered 0/i, increases, matching the trend that lowering the wall temperature increases the liquid impingement 
effect. The physical reasoning for this is that as the wall temperature is lowered, there is less production of vapor at 
the surface to resist the momentum of the liquid droplets approaching the wall. At lower values of T w , more droplets 
are able to penetrate the vapor film and make contact with the tube surface, effectively increasing the heat flux 
beyond normal film boiling values. 9fb is a maximum value of 1 at T w = T wet . For wall temperatures below T we , the 
correlation is no longer used because the flow has entered transition boiling. It was estimated that a power of 3 for 
Ofl, was appropriate. 


B. Nucleate Boiling HTC Correlation 

Many correlations in the literature for flow boiling in the nucleate boiling regime utilize the boiling (Bo) number 
Bo = q /(Gy/ V ). 13-17 These correlations are based on steady-state heated tube data with the main application of 
electronics cooling or controlled heating where the heat flux is applied and known. In these cases the variable of 
interest is the wall superheat, i.e. how hot the tube wall material will get for a known heat flux and known mass flux 
of the liquid. In chilldown, the heat flux is the variable of interest, or the unknown variable. The wall superheat is 
the known variable since it is known at each time step of the simulation. If one were to use one of the many 
correlations in the literature with the Bo number as one of the parameters then this would require iteration until the 
heat flux on the LFIS of the equation matched the heat flux on the RHS. This has the potential of significantly 
slowing down the simulation run time. Additionally, it has been shown that the transient nature of quenching (i.e. 
chilldown) data causes it to deviate from steady-state data. 18 

Therefore the approach taken in the current work is to find a correlation that is a function of all the variables that 
are known, e.g. the mass flux, the fluid and wall properties, etc. The correlation developed by Chen 19 provides the 
base form for a very robust correlation. It is based on the physics of the heat transfer process and thus does not 
involve an empirical correlation in terms of the Bo number. The Chen correlation will be described here. 

The nucleate boiling regime HTC, h„b, is the sum of the HTC responsible for forced convection with the two- 
phase fluid, hf c , and the HTC responsible for the heat transfer that arises from bubble nucleation and growth in 
surface cavities, hi,. The former is sometimes called macroconvection whereas the latter is called microconvection. 
The equation for h„b is: 


K b = h fc+K 


( 2 ) 


where the two components are: 


/ Vf =0.023[^(l-xJ]°V 4 |^ 


h, = 0.00122 


i 0.19 si0A5 0.49 

K l C p,l Pi 

0.5 0.29 0.24 0.24 

Pi Ylv Pv 


[K-TjP)f 24 [P at {T w )-pf 75 S 


( 3 ) 

( 4 ) 
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All fluid properties are calculated at the saturation temperature. Eq. (3) is just the single-phase liquid HTC from 
the Dittus-Boelter equation with the factor F to account for the difference between the liquid Reynolds (Re) number 
and the actual two-phase Re number. The difference between the two-phase thermal conductivity and the liquid 
thermal conductivity are considered negligible since k, » K. Also the two-phase Prandtl (Pr) number is taken 
simply as the liquid Pr number since both Pr / and Pr v are typically very close for most fluids. F was found 
empirically by Chen as: 


F = 


1 


\ 


— + 0.213 


( 5 ) 


The empirical term S in Eq. (4) accounts for the difference in the actual superheat and total superheat seen by 
bubbles growing from surface cavities. The difference is caused by the presence of the thermal boundary layer. S 
was found by Chen to be: 


S = 


1 


( 6 ) 


l + 2.53xl0^7?e / U7 F L4625 

Eqs. (2)-(6) represent the correlation used by this current work to determine the HTC in the nucleate boiling regime 


C. Transition Boiling HTC Correlation 

Transition boiling consists of a mix of film and nucleate boiling, it was thought that the correlation could use 
either or both of the film boiling and nucleate boiling correlations with an additional empirical constant. Since 
nucleate boiling is an order or two magnitude greater than film boiling, a fit was applied to the data with the Chen 
correlation HTC and the nondimensional temperature, 9,b = ( T wet - 7j,)/( T wet - T sal ). The final equation for the 
transition boiling HTC was: 

h tb =Q.52Wr\ b ( 7 ) 

The nondimensional temperature 9,b is always between 0 and 1. Hence, in transition boiling the heat transfer is only 
a fraction of the full nucleate boiling heat transfer because there is some presence of film boiling. When 7), is larger 
than T wet there is no nucleate boiling yet and only film boiling. As T w decreases below T wel , 9,b becomes nonzero and 
grows larger, indicating more nucleate boiling heat transfer is present. 

D. CHF Correlation 

The CHF has been reported to be a strong function of the We number based on the length downstream of the 
inlet. 20 The following correlation that agreed well with the data was used in the current work: 

q CHF = 0.0527 Gy h We^ (8) 


where We. = G 2 zl(pia). 

E. Rewetting Temperature Correlation 

The model developed by DeSalve and Panel kr is used to predict the rewetting temperature. It calculates the 
rewetting temperature as: 


T wet = T al (P) + 0.29d l [T ms -T m (/>)](l + 0.279G a49 ) (9) 

where d = exp[3.06xl0 6 /(£,/i s C' AS )]erfc[1751.5(A; J /i i C ; , s ) 0 ' 5 ] takes into account the effect of the contact temperature 
between the liquid and the solid. T us is the theoretical maximum superheat calculated from Spiegler et al. 22 : 

T ms = 0.844T r (10) 
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where T cr is the critical temperature of the fluid in Kelvin. Both models are also described in Carbajo. 2 ’ 

F. ONB Temperature Correlation 

A simple empirical fit to the data was used to correlate the ONB temperature: 

T ma =T Sttt {P) + 0.001\P + 5 (11) 

Again, the temperature is in Kelvin. This correlation is more empirical and likely less applicable to other fluids 
than one would hope. Therefore future work for this group will include developing a much more robust correlation 
for Tqnb • Eq. (11) will suffice for this paper since the main goal here is to show the importance of the accuracy of 
the other more robust correlations for T wel , hjb, h„b, h,b, and the CHF. 

VI. ID Simulation Description 

The simulation of the chilldown of the test section tube was developed in MATLAB . The test section, 
(geometry shown in Figure 2) was divided into lumped-parameter nodes and the finite volume technique was used to 
solve the transient energy equation in the axial direction. It was assumed that the temperature gradients in the radial 
and azimuthal directions were negligible in comparison to the gradients in the axial direction. This is appropriate for 
thin tubes like the one used in the current work. Under this assumption, the energy equation for the solid tube 
becomes: 


ST (y T 

Ps C p,s -^ =k s -JgT + 9am, ( Z ) + q'para.mc ( Z ) ( 12 ) 

where the subscript “s” refers to a property of the solid tube. There were two external heat fluxes: one from the fluid 
convection, q conv , and another from parasitic heat, q parasitic- Radiation and gas conduction from the surroundings 
outside the test section tube made up the parasitic heat and were estimated for the simulation. A central difference 
scheme was used for the second-order spatial derivative and the first-order fully-implicit Euler method was used for 
time integration. Using m = pA r ,Az. q con v = h conv (Tf- 7’,,), and the fact that the convective heat flux acted over the 
inside surface of the tube and the parasitic heat flux acted over the outside surface of the tube, Eq. (12) for each node 
i becomes 


mC ( T j+1 -T j .) = (kAAt/Az)( T j t. - 2 T J+1 + T J+ \ ) 

p,s y w,i w,i J \ s cs J y vty+1 w,i w,i — 1 j 


+ /rDAzAth . 


■ ( T u ~ ) + n{D + 2t tuhe ) AzAt (13) 


Where j and j+1 represent the old and new times, At is the time step, A z is the node length, m is the mass of each 
node, t tu b e is the thickness of the tube, A cs = 7t/4[(D + 2t lube ) 2 - D 2 ] is the cross-sectional area of the tube, and 7} is the 
fluid temperature. Figure 3 illustrates Eq. 10 for node i. Note that each solid node consists of the entire cross section 
and therefore envelopes the interior fluid node. 

The parasitic heat flux was a summation of the heat flux due to the conduction between the thin gas in the 
vacuum chamber and the tube outer surface, q gc , and the heat flux due to radiation between the vacuum chamber 
inner surface and the tube outer surface, q ra d- For two concentric cylinders with a gap in between filled with 
motionless air, the air having thermal conductivity of K e , the equation for q becomes: 

^=^(^-7;) cosh- 1 


(D v J+jD + t tu J 
D vac (D + t, uhe ) 


D + 1„ 


(14) 


where T wv ac is the temperature of the vacuum chamber taken as a constant 293 K value, and D vac = 6.02 cm is the 
inner diameter of the vacuum chamber. The thermal conductivity of air is a strong function of the pressure and is 
estimated by: 24 
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K _ K , CT a mb I 

K, 1 ^[^-(a+Oll 


( 15 ) 


where K 0 is the thermal conductivity of air at atmospheric pressure which is approximately 26.3 mW/mK, P vac is the 
vacuum pressure which was held at 14 Pa, T amb is the ambient temperature taken as 293 K, and c is a constant equal 
to 7.6 X 10' 5 with units of Nm 'K' 1 . The q rac i portion is estimated by an equation describing the radiation heat flux 
between two infinitely long concentric cylinders: 


tfrad 




(1 ^w.vac ) 


w.tube 


D + tfube 
V D vac J 


(16) 


Vacuum 


where a SB is the Stefan-Boltzmann constant, c wtuhe is the emissivity of the tube, and e w vac is the emissivity of the 
vacuum chamber. The emissivity for each surface was estimated as 0.45 since both are stainless steel. 

Eq. 10 was solved simultaneously for all nodes by forming a tridiagonal matrix and solving by the Thomas 
algorithm. For each test, the test section was discretized into 40 nodes. With a test section length of 57.2 cm, this 
gave a node spacing of A z = 1.34 cm. This node spacing and a time step of At = 0.01 s were found to be small 
enough to achieve convergence. The boundary conditions were as follows. The first node, located at the test section 
inlet, was assumed to start and remain at the inlet liquid temperature, determined from a TC placed just upstream of 
the bypass valve (see Figure 2). The last node was assumed to have the same heat flux as the node prior to it. It 
should be noted that different boundary conditions than these were attempted but did not significantly change the 
results of the simulation. This is because the convective heat flux with the fluid dominated over the conduction 
between nodes for the entire test except after the ONB temperature. Also, the last node was several nodes 
downstream of the locations of interest, which are TC-1 and TC-2 in Figure 2. Initially, all nodes except the first 
node were set to the first temperature reading of the given test, which was typically around 293 K. 

This simulation modeled the energy equation 
of the solid. The model assumed a constant mass 
flow rate and time-independent pressure 
distribution of the fluid during the chilldown 
process, and therefore the momentum and 
continuity equations for the fluid did not need to 
be solved. In reality there were some oscillations 
in mass flow rate and pressure for each test. An 
average value of the mass flow rate data was used 
to simulate a given test. The pressure drop through 
the test section was assumed to be linear. The two 
time -averaged pressure readings from PT-1 and 
PT-2 shown in Figure 2 were used to calculate a 
linear pressure drop along the test section. The 
average mass flow rate and fluid pressure found 
from each test were used in evaluating h conv . The 
fluid temperature and quality, however, had to be 
found at each time step and each node. The 
method used can be described by referring to 
Figure 3. Next to each solid node there is a 
contiguous fluid node. At each time step the 
enthalpy of a fluid node, y„ was calculated by: 
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Figure 3. Numerical Simulation Setup. 


Ti = Yt - 1 + Wi [ 4 1 (^ 2 g)] 


(17) 


where y bl and q CO nvj-i are the enthalpy and the heat input in the previous node. For the first node, the enthalpy was 
defined by the inlet temperature and inlet pressure measured just upstream of the bypass valve. With the enthalpy 
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Table 1 - L-norm of the 
temperature errors 


defined at each fluid node, the equilibrium quality was found by x e = (y t - yi_ sat )/yi v where yi iSa , and y iv are the saturated 
liquid enthalpy and the heat of vaporization at the local pressure. The temperature was also found knowing the 
enthalpy and the pressure for each node. With the local temperature, quality, and pressure, all of the fluid properties 
were determined at each node using Refprop, which enabled the use of the HTC correlations to evaluate h com , at each 
node. 

Even though an implicit formulation was used in Eq. (13), the enthalpy of 
each fluid node is solved explicitly in Eq. (17). Therefore a convergence study 
was performed to show that the chosen time step was small enough to achieve 
convergence. Time step values of At = 0.04 s, 0.02 s, 0.01 s, and 0.005s were 
used to simulate the fastest chilldown run of all 55 tests. The predicted 
temperature traces at TCI are shown in Figure 4 for each time step. It is difficult 
to distinguish the difference between the chilldown 
curves, indicating an acceptable level of convergence. 

For further measure, the L 2 norm of the errors for At = 

0.04 s, 0.02 s, and 0.01 s were found by assuming the At 
= 0.005 s solution was exact. The results are shown in 
Table 2. When changing from At = 0.02 s to A? = 0.01 s 
the L 2 norm decreases by a much larger fraction than 
when changing from A? = 0.04 s to At = 0.02 s. Therefore 
a value of At = 0.01 s was used for all simulations. 

The logic to determine what correlation to use - 
either film boiling, transition boiling, or nucleate boiling 
- went as follows: 


Time step, s 

L 2 -norm 

0.04 

31.6 

0.02 

19.0 

0.01 

7.4 


• Calculate T wet for the mass flow rate and local 
pressure 

• If T w ^ T wet then h conv kfb 

• If T w < T wet , then: 

o If T w > T onb , then: 

- Calculate CHF 

■ Calculate q nh prediction from 

Chen correlation 
" If q^nb ''' 7 / v ib then h conv h tB 

" i f 7 nh ^ 7 chf then h conv h nB 

• If T w < T onb then h conv = h DB where h DB is the single phase Dittus-Boelter HTC 



Figure 4. Predicted chilldown curves of simulation for 
several time step values. 


VII. Simulation Results 

The simulation was run for 55 different tests with ranges of 
mass fluxes from 61.2 to 1150 kg/m 2 s, inlet pressures from 175 to 
817 kPa, and subcooled inlet temperatures from 0 to 14 K below 
the saturation temperature. Plots of the measured and predicted 
chilldown curves at the upstream TC stations for six different mass 
fluxes are shown in Figure 6-Figure 8. The mass flux, inlet liquid 
Re number, and inlet pressure are listed in the caption of each plot. 

As can be seen from these plots, the simulation performed well over all conditions. The main cause of the slight 
disagreement at all mass fluxes is likely due the simplification of using an average mass flow rate and a time- 
independent pressure distribution when there were actually some fluctuations in mass flux and pressure 
measurements for each test. An appropriate measure of the usefulness of the correlations of the current work should 
use the chillldown time, or how long it takes for the wall to reach the temperature of the fluid, since this is one of the 
main concerns from an engineering standpoint. Some of the tests were not able to reach the saturation temperature 
because at such low temperatures the parasitic heat began to overcome the convective heat flux, so two different 
time errors were considered - 1) the time from start to the rewetting temperature, t wet , and 2) the time from start to 
the ONB, t 0NB . Both the mean absolute error (MAE) and the mean absolute percent error (MAPE) for t wet and t ONB 
are shown in Table 2. The two errors are defined as MAE = 1 IN x 2(|? Jlm - t data |) and MAPE = 1 IN x £(|? Jim - 


Table 2 - Percent MAE of temperature 
prediction in each boiling regime 


Time Metric 

MAE, s 

MAPE, % 

twet 

4.83 

24 

t()NB 

5.15 

25 
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tdata! tdata) x 100%, where the subscripts t sim and t data stand for the time predicted by the simulation and the time 
measured from the data. The summation is over the N total number of data points in each regime for all 55 tests at 
both TC stations. Times t sim and t data are predicted to within 25%. 
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VIII. Conclusion 

The numerical simulation of the chilldown of a short, thin, vertical tube with LN 2 flowing downward in the 
current work was able to predict the transient temperature of the tube with good agreement with the data over a wide 
range of conditions. These results support the validity of the HTC correlations proposed in the current work for all 
three of the boiling regimes (even for nonequilibrium, i.e. negative quality values, at the inlet), the correlations for 
the dividing points of the boiling regimes (the rewetting temperature, the CHF, and the ONB temperature), as well 
as the logic employed when choosing which correlation to use under local conditions. These correlations could be 
used to improve the accuracy of commercial simulation packages. Future work includes: a full description of the 
experimental apparatus and the experimental results which will include multiple flow angles in addition to the 
downward flow angle considered in the current work, further analysis and data comparison of the correlations in the 
current work, and extension of the correlations in the current work to different flow angles. 
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